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ABSTRACT 

Our  concern  is  to  "unfold"  and  flatten  the  curved,  convoluted  sur- 
faces of  the  brain  in  order  to  study  the  functional  architectures  and 
neural  maps  embedded  in  them.  In  order  to  do  this,  it  is  necessar>'  to 
solve  the  general  map  makers  problem  for  representing  curved  surfaces 
by  quasi-isometric  planar  models.  This  algorithm  has  applications  m 
areas  other  than  computer  aided  neuroanatomy,  such  as  robotics  motion 
planning  and  geophysics. 

The  algoritlim  we  have  written  maximizes  the  goodness  of  fit  of  dis- 
tances in  these  surfaces,  to  those  in  a  planar  configuration  of  points. 
We  illustrate  this  algorithm  with  a  flattening  of  monkey  visual  cortex, 
which  is  an  extremely  complex,  folded  surface.  We  find  distance  errors 
in  the  range  of  several  percent,  with  isolated  regions  of  larger  error,  for 
the  class  of  cortical  surfaces  which  we  have  so  far  studied. 
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INTRODUCTION 

The  map-maker's  problem  is  to  find  a  flat  representation  of  a  curved  surface:  for 
example,  the  surface  of  the  earth.  Classical  map  making  has  been  restricted  to  the 
relatively  simple  spherical  surface  of  the  earth.  In  the  case  that  the  surface  of  interest 
is  complex,  and  possibly  non-convex,  there  are  no  known  methods  of  finding  quasi- 
isometric  planar  representations. 

The  solution  to  this  problem  is  of  importance  to  computer  aided  neuro-anatomy, 
since  it  is  often  desired  to  view  the  surface  of  various  cortical  areas  in  a  planar  model. 
Primate  cortex  is  highly  convoluted,  and  provides  one  of  the  more  complex  surfaces 
encountered  in  practical  applications. 
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Motion  planning  in  robotics  is  another  area  of  application  for  which  the  finding  of 
shortest  distances  on  polyhedral  surfaces  is  of  importance  (Sharir  and  Schorr,  1984). 
Other  areas  of  bio-physics  and  geo-physics  would  seem  to  provide  possible  areas  of 
application  of  a  generalized  map  makers  algorithm. 

Our  interest  is  in  obtaining  a  flat  representation  of  the  cortical  surfaces  of  the 
brain,  because  the  detailed  maps  of  sensory  and  other  neural  data  embedded  in  these 
surfaces  are  easiest  to  study  and  measure  when  they  are  shown  on  a  plane  surface. 

This  report  describes  the  method  we  use  to  find  an  optimal  quasi-isometry  that 
maps  an  arbitrary  curved  surface  into  a  plane.  This  mapping  is  optimal  in  the  sense 
that  it  is  derived  from  a  variational  principle  that  optimizes  the  overall  fit  between  dis- 
tances in  both  the  curved  and  planar  surfaces.  The  mapping  is  a  quasi-isometry 
because  it  optimizes  the  fit  of  distances,  over  multiple  scales,  rather  than,  say,  local 
angles  (in  which  case  it  would  be  quasi-isothermal,  or  conformal). 

We  use  a  quasi-Newton-Raphson  minimization,  applying  it  to  the  distance  matrix 
of  inter-point  distances  up  to  some  order  of  distance  from  each  node.  Minimal  geo- 
desic distances  are  calculated  in  the  surface,  and  these  form  the  initial  distance  matrix. 
The  planar-distance  matrix  is  initially  generated  from  random  numbers.  Planar  points 
are  moved  by  the  gradient  of  a   "stress,"  which  measures  the  quality  of  the  isometry. 

We  have  found  that  different  random  starting  configurations  converge  to  virtually 
identical  final  states.  This  finding  indicates  lack  of  "trapping"  in  local  minima,  and 
also  indicates  the  robustness  of  our  method.  We  also  show  that  in  general,  local  (i.e., 
nearest-neighbor-only)  distances  are  insufficient  for  this  problem.  We  base  this  con- 
clusion on  general  arguments  and  on  our  own  experience  with  local-distance  problems. 

Description  of  the  algorithm 

Given  an  arbitrary  cur\'ed  surface,  we  wish  to  "flatten"  it.  Because  surfaces  gen- 
erally have  non-zero  Gaussian  curvature,  there  is  usually  no  isometrN'  between  the  sur- 


face  and  the  corresponding  plane  (do  Carmo,  1975).  In  other  work  (Kaplow  and 
Schwartz,  1986),  we  describe  ways  to  measure  the  mean  and  Gaussian  curvature  of  a 
polyhedral  surface. 

We  state  the  map-maker's  problem  in  the  following  way:  Given  a  polyhedral 
surface,  determine  the  entire  set  of  interpoint  distances.  This  lower-triangular  distance 
matrix  is  of  size  N(N-l)/2  for  N  nodes.  We  then  determine  a  set  of  N  points  in  the 
plane  such  that  the  corresponding  distance  matrix  in  the  plane  provides  a  best  fit,  in 
the  least-squares  sense,  with  the  original  distance  matrix  of  the  surface.  This  pro- 
cedure should  yield  an  optimal  quasi-isometric  mapping  of  the  surface  into  the  plane. 

Distances  in  polyhedral  surfaces 

One  difficulty  with  this  algorithm  is  that  until  recently,  no  algorithms  were  known 
that  can  find  distances  on  arbitrary'  (i.e.,  non-convex)  polyhedral  surfaces.  Sharir  and 
Schorr  (1984)  describe  an  algorithm  that  finds  minimal  distances  in  a  convex 
polyhedral  space,  and  other  researchers  (Mount,  1984,  Mitchell  and  Papadimitrious, 
1984, 0'Rourke  et  al,  1984)  have  described  algorithms  to  find  distances  on  non-convex 
polyhedrons.  These  algorithms  are  complex,  and  we  know  of  no  actual  implementa- 
tions of  them.  Elsewhere  we  describe  an  algorithm  we  have  developed  that  finds 
minimal  distances  in  arbitrary  polyhedral  surfaces.  In  the  present  paper,  we  assume 
that  we  have  obtained  the  distance  matrix  of  the  original  surface  by  applying  that  algo- 
rithm to  compute  the  discrete  minimal  distances  (Wolfson  and  Schwartz,  1987). 

Variational  algorithm 

Our  mapping  problem  is  ver)'  similar  to  the  conventional  "multi-dimensional" 
scaling  problem,  which  is  used  to  optimally  map  point  sets  from  N  to  M  dimensions, 
in  which  M  is  usually  1,  2  or  3,  and  in  which  N  is  usually  large  (e.g.,  10-50).  This  is 
a  common  cluster-analysis  procedure  that  is  often  used  to  visualize  multivariate  data 
sets  (Sammon,  1969,Schiffman  et  al.,  1981). 


In  the  present  case,  because  we  are  mapping  between  two-dimensional  surfaces, 
we  would  expect  good  performance,  especially  because  our  surfaces,  although  highly 
convoluted  or  folded,  do  not  usually  represent  "crumpled"  spheres,  but  rather  crumpled 
sections  of  spheres.  In  other  words,  our  surfaces  are  not  closed;  and,  as  indicated  by 
separate  measurements,  their  integral  mean  curvature  is  not  very  large  (Kaplow  and 
Schwartz,  1986). 

Description  of  the  algorithm 

Let  the  (minimal  geodesic)  distance  in  the  surface  be  organized  as  a  lower- 
triangular  matrix  d^  Let  the  (unknown)  distances  in  the  plane  be  organized  as  a 
lower-triangular  matrix  J,..  We  pose  the  map-maker's  problem  as  minimizing  the 
least-squares  goodness  of  fit  between  these  two  distance  matrices.  Thus,  we  propose 
to  minimize  the  quantity  L  below: 


1  i=N  kj  -  ^ijj 

L  =  -^£-5-^ ^  (1.0) 

c  i<j         dij 


c=idij  (1.1) 

Because  the  distance  matrix  contains  all  possible  interpoint  distances  in  this 
discrete  problem,  an  optimal  preservation  of  the  distance  matrix  is  equivalent  to  an 
optimal  preservation  of  the  metric  structure  of  the  original  surface,  in  a  planar 
representation. 

We  use  the  gradient-descent  variational  technique.  Starting  with  a  random  planar 
point  configuration,  we  calculate  the  initial  distance  matrix  from  these  random  points. 
Then  we  analytically  calculate  the  gradient  of  equation  (1.0)  above,  and  move  in  the 
direction  of  the  descending  gradient  by  an  amount  scaled  by  the  Hessian  of  equation 


(1.0).  This  is  the  quasi-Newton  Raphson  technique.  The  quantity  L  in  equation  (1.0), 
which  represents  the  goodness  of  fit  of  the  metric  structure  of  the  data  and  planar-point 
sets,  we  term  the  "stress."  Our  approach  minimizes  this  stress.  Plots  of  stress  as  a 
function  of  the  number  of  interations  give  some  indication  of  how  well  this  algorithm 
performs.  In  addition,  testing  the  results  with  several  random  starting  configurations, 
as  shown  in  figure  2,  can  be  used  if  there  is  any  doubt  about  the  validity  of  the  solu- 
tion. 

Experiments  with  neighbor  and  next-neighbor  distances 

We  implemented  the  algorithm  of  equation  1.0  above  before  we  had  finished  its 
companion  geodesic-distance-determination  algorithm  (Wolfson  and  Schwartz,  1987). 
Thus,  our  early  experiments  involved  only  neighbor  distances  (i.e.,  the  edge-lengths  of 
the  original  polyhedron)  and  next-neighbor  distances  (which  we  obtained  easily  by 
applying  the  law  of  cosines  across  two  neighboring  triangles).  The  results  of  these 
experiments  were  disappointing.  When  the  algorithm  was  initialized  with  a  random 
starting  configuration,  there  generally  was  no  convergence  to  the  expected  final  state. 
Instead,  the  configurations  we  obtained  seemed  to  be  "folded."  In  other  words,  large 
regions  had  the  right  structure  but  were  in  the  wrong  position  in  the  final  planar 
configuration. 

This  result  was  easy  to  understand.  When  only  short-range  distances  are  avail- 
able, little  penalty  is  imposed  on  globally  incorrect  patches  whose  internal  structural 
details  are  correct.  This  is  because  the  penalty  is  essentially  imposed  only  on  the 
boundaries  of  the  regions.  The  boundaries  of  patches  of  "diameter"  m  have  a  weight 
of  0(m),  whereas  the  locally  correct  interior  structure  has  a  weight  of  O(m-). 

The  only  condition  under  which  we  were  able  to  obtain  successive  flattenings 
with  short-range  distances  was  by  supplying  the  algorithm  with  an  initial  state  that  had 
been  relatively  well  flattened  beforehand  (i.e.,  through  human  intervention).  Obvi- 
ously, such  a  procedure  cannot  be  called  a  flattening  algorithm.    However,  by  adding 


longer-range  distance  terms,  we  successfully  flattened  polyhedra  from  random  starting 
configurations,  as  described  below. 

Heuristic  simplication  of  the  algorithm 

The  size  of  the  distance  matrix  (and  hence  the  complexity  of  this  algorithm) 
scales  as  N",  where  N  is  the  number  of  nodes.  The  computer  memory  and  CPU  time 
required  for  large  N  are  therefore  prohibitive.  We  restricted  the  size  of  the  neighbor- 
hood of  distances  represented  for  each  node.  This  restriction  provided  a  heuristic  in 
which  a  "patch"  of  surface  around  each  node  contributed  to  the  problem.  Thus,  our 
distance  matrix  became  relatively  sparse.  By  storing  this  data  in  a  forest  of  binar>- 
search  trees,  rather  than  a  matrix,  we  achieved  reasonable  performance  on  a  Sun  2 
microprocessor  system  (which  is  roughly  comparable  to  a  VAX  750).  The  extent  of 
these  patches  is  an  empirical  problem.  In  our  experience,  a  patch  whose  size  is 
roughly  10%  of  the  surface  area  of  the  entire  problem  yields  good  results.  Thus,  for 
example,  as  shown  in  Figure  3d,  with  0(1000)  nodes,  a  patch  of  "diameter"  of  about 
10  gave  good  results. 

Examples  of  surface  flattening 

Example  1:   Hemisphere 

We  generated  a  hemispherical  shell,  calculated  the  distance  matrix  analytically 
(for  an  eight  neighborhood),  and  "flattened"  it  with  our  algorithm.  The  result  is  shown 
in  Figure  1. 

Example  2:    Visual  cortex 

The  surface  of  the  visual  cortex  of  a  macaque  monkey  was  digitized,  triangulated. 
and  flattened  as  in  Figure  2.  Figure  2  shows  the  "opercular"  surface  of  visual  conex. 
The  operculum  (Latin  for  "roof)  is  an  easily  accessible  part  of  the  visual  cortex.  In 
fact,  as  shown  in  Figure  3b,  the  opercular  cortex  can  be  flattened  physically  with  rela- 
tively  little   distortion.     We   studied  the  error  term   from  physical   flattening  of  the 


opercular  surface  of  cortex  by  making  India-ink  dots  on  cortex  which  we  then  flattened 
and  filmed  with  a  movie  camera.  Figure  3b  is  an  picture  from  this  series.  For  the 
opercular  surface,  once  the  isotropic  expansion  has  been  corrected,  the  physical  flatten- 
ing shear  error  is  typically  5%  to  10%.  Our  digital  flattening  errors(  for  the  opercu- 
lum), on  the  other  hand,  are  typically  in  the  1%  range. 

The  usefulness  of  digital  flattening  is  clearer  when  we  consider  the  entire  striate 
cortex  instead  of  just  the  operculum.  Figure  3c  is  a  computer-generated  reconstruction 
of  the  entire  cortex,  viewed  from  below.  Although  the  cortex  is  quite  convoluted,  this 
view  shows  that  its  surfaces  are  generally  rather  flat.  In  fact,  it  is  easy  to  fold  a  sheet 
of  paper  into  the  petal-like  shape  shown  in  Figure  3c  without  tearing  or  stretching  the 
paper. 

Figure  3d  shows  the  flattening  of  the  entire  cortex.  The  model  in  this  example 
was  a  polyhedron  consisting  of  about  2500  triangles.  We  flattened  the  model  by  using 
a  12-neighbor-deep  distance  around  each  node.  Typical  errors  (  defining  error  as  the 
average  local  difference  of  edges  in  3  dimensions  and  2  dimensions)  was  typically  in 
the  range  of  a  few  percent. 

Other  applications 

As  noted  earlier,  this  algorithm  is  a  solution  to  the  general  map-maker's  problem. 
The  visual  cortex  is  the  largest  and  arguably  the  most  convoluted  of  the  cortical 
regions,  and  our  cortical  data  seems  to  be  rather  more  complex  than  the  data  in  other 
common  instances  of  the  problem.  The  good  results  we  have  obtained  with  this  data 
encourage  us  to  think  that  our  methods  can  deal  competently  with  a  variety  of  prob- 
lems related  to  the  mappmg  of  surfaces.  In  particular,  we  believe  that  problems  in 
robot-motion  planning  and  autonomous  vehicle-navigation  might  yield  to  an  applica- 
tion of  this  alcorithm. 
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Figure  Captions 

Figure  1.  A  hemispherical  shell  was  flattened  by  our  algorithm.  This  figure 
shows  the  3D  wire-frame  model  of  the  original  surface  and  also  the  flattened  version 
of  the  surface.  A  "neighborhood"  extending  at  least  eight  nodes  from  each  node  deter- 
mined the  "patch"  of  surface  for  which  we  computed  the  distances. 

Figure  2. 
A  section  of  the  "operculum"  of  visual  cortex,  reconstructed  from  serial  sections  of 
monkey  brain,  was  flattened  by  our  algorithm.  The  three  different  versions  of  the 
planar  surface  were  obtained  in  successive  runs,  with  different  random  starting 
configurations.  Because  the  algoritlim  does  not  constrain  the  angle  of  rotation,  the 
solution  is  determined  only  up  ro  a  rotation  and  reflection.  However,  size  is  preserved. 
The  three  solutions  are  vir^-ally  identical,  suggesting  that  no  local  minima  that  might 
be  present  have  obscured  the  solution.  This  suggestion  is  reinforced  by  the  plot  of  the 
stress  as  a  function  of  iterations,  shown  at  the  lower  right.  In  this  parameter  (which 
determines  the  overall  least-squares  goodness  of  fit)  there  is  some  large  oscillation 
early  in  the  run;  but  after  about  100  iterations,  all  three  runs  converge  to  the  same 
stress  and  to  the  same  geometric  solution.  The  curved  appearance  of  these  figures  is 
an  illusion  caused  by  the  curved  nature  of  the  original  polyhedron.  In  fact,  all  three  of 
these  examples  are  of  course  planar  representations. 

For  this  model,  minimal  distances  were  calculated  usuing  the  algorithm  of  Wolf- 
son  and  Schwartz  (1986).  Use  of  an  eight-node  neighborhood  of  the  600-node  model 
limited  the  size  of  the  distance  matrixes. 

Figure  3a. 
NOTE;  ORIGINAL  FIGURE  IS  GRAY  SCALE  FOR  JOURNAL  REPRODUCTION. 
A  block  of  monkey  brain.    The  smooth  roof  (operculum)  of  visual  cortex  is  the  bullet- 
shaped  region  at  the  left. 
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Figure  3b. 
NOTE:  ORIGINAL  FIGURE  IS  GRAY  SCALE  FOR  JOURNAL  REPRODUCTION. 

The  opercular  visual  cortex  is  physically  flattened  by  being  pressed  between  glass 
plates  while  being  filmed  with  an  8mm  movie  camera.  The  India-ink  dots  visible  on 
the  cortex  were  used  to  measure  the  distortion  caused  by  physical  flattening. 

Figure  3c. 
NOTE:  ORIGINAL  FIGURE  IS  GRAY  SCALE  FOR  JOURNAL  REPRODUCTION. 

A  computer-generated  reconstruction  of  the  entire  striate  cortex,  viewed  from  the  bot- 
tom. The  opercular  surface  is  away  from  the  viewer,  while  the  convoluted  "calcarine" 
cortex  appears  as  a  flower-petal  shape  closer  to  the  viewer. 

Figure  3d.  The  entire  striate  cortex  from  Figure  3c,  as  flattened  by  our  algorithm. 
Greater  densities  of  triangles  in  the  model  indicate  the  regions  of  higher  curvature  on 
the  cortical  surface. 
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